A Quick Simulation Method for Aero-Optical Effects Based on a Density Proxy Model

Aero-optical effects caused by high-speed flow fields will interfere with the transmission of starlight, reduce the accuracy of optical sensors, and affect the application of celestial navigation on hypersonic vehicles. At present, the research of aero-optical effects relies heavily on the flow field simulation of computational fluid dynamics (CFD), which requires a great deal of computing resources and time, and cannot satisfy the demand of the rapid analysis of aero-optical effects in the engineering design stage. Therefore, a quick simulation method for aero-optical effects based on a density proxy model (DP-AOQS) is proposed in this paper. A proxy model of the turbulent density field is designed to replace the density field in the CFD simulation, and the proxy model is parametrically calibrated to simulate the optical characteristics of the turbulent boundary layer (TBL) in the external flow field of the optical window. The performance of DP-AOQS in the visible light band is verified from the perspectives of density field distribution, optical path difference (OPD), and fuzzy star map. The simulation results show that the method can quickly provide the distortion results of aero-optical effects in different flight conditions on the premise of ensuring the simulation accuracy. The research in this paper provides a new analytical method for the study of aero-optical effects.


Introduction
With the rapid development of the aerospace industry in various countries, hypersonic vehicles have become a major research focus of aircraft technology [1,2]. In order to improve the autonomy and accuracy of navigation and guidance, hypersonic vehicles are often equipped with optical imaging detection equipment to provide detection information for navigation and guidance. The application of star sensors in celestial navigation is a particularly extensive research direction at present [3,4]. However, when hypersonic vehicles fly at a high speed in the atmosphere, complex flow structures, such as turbulent boundary layers, shock waves, and expansion waves and mixed layers, may be formed outside the optical window, which makes the density field, temperature field, and velocity field of the air change dramatically [5,6]. When light passes through the complex flow field, the high-speed flow field will disturb the transmission of light and cause optical distortion; that is, aero-optical effects which can make the target image formed by imaging detection equipment appear subject to energy attenuation, offset, blur, jitter, and other phenomena [7][8][9][10]. Therefore, the research of aero-optical effects is of great significance for the application of optical sensors in hypersonic vehicles.
The law of light disturbance caused by high-speed flow field is important to lay a foundation for schemes designed to reduce the influence of aero-optical effects. At present, the main method of research is to obtain the external flow field data of hypersonic vehicles by using high-precision numerical calculation methods, such as computational fluid dynamics (CFD) and direct simulation Monte Carlo (DSMC) [11][12][13]. Then, the optical distortion and image quality evaluation of optical imaging equipment can be obtained through geometric optics, physical optics, statistical optics, and other optical transmission calculation methods [14][15][16]. The focus of those methods is mainly on how to obtain the flow field data, which requires a huge amount of computation and time. When analyzing the influence of flight conditions on the aero-optical effects, most of them adopt the method of repeatedly solving the flow field to calculate the aero-optical effects under the corresponding conditions [17][18][19]. Some scholars have studied the vortex structure of the boundary layer flow field by the CFD method, but the solution of complex differential equation consumes a lot of computational resources [20]. Although the CFD and DSMC methods have a high accuracy in calculating the density field, and the aero-optical effect results calculated from them have a high reliability, the calculation takes a long time and the utilization rate of the data is low. Due to the complexity of turbulence, it is necessary to recalculate the full flow field data when the flight conditions change slightly. Numerical simulation methods such as these, with a high precision and high computation, are used to calculate the density field to analyze the aero-optical effects, which cannot satisfy the requirement of rapidity when conducting aero-optical engineering research at the aircraft design stage. Although some scholars have improved various numerical methods based on hydrodynamics [21], in terms of the aero-optical effects, we mainly focus on the optical properties of the flow field. Therefore, it is necessary to find a quick simulation method for the aero-optical effects expected for the research method based on CFD.
As the calculation of a complex flow field is the main factor affecting computation in aero-optical simulations, a proxy model of the flow field can effectively improve the simulation efficiency of aero-optical effects. Many scholars have verified, through many experiments and numerical simulations, that the aero-optical effects have a corresponding relationship with the boundary layer thickness, Mach number, incoming flow density and other parameters of the optical window boundary layer, and obtained scaling law models such as the root mean square of the optical path difference and jitter angle [22][23][24][25][26]. These high-fidelity aero-optical distortion laws are the basis for matching the flow-field proxy model with the actual high-speed flow field.
The contribution of this paper is to design a boundary layer density proxy model to replace the density field in CFD numerical simulation and realize the rapid simulation of aero-optical effects, as shown in Figure 1. In this paper, the density proxy model is based on the density ellipsoidal vortex model in the aero-optical effects steady-state simulator proposed previously [27]. In contrast to the previous method, this paper considers the density distribution inside the spherical vortex and converts the large-scale vortex structure into the superposition of multiple ellipsoidal vortices. To conform to the real high-speed flow field environment, the motion model of the ellipsoidal vortex structure is added on the basis of the previous static structure to make the density proxy model closer to the real high-speed flow field. In order to ensure the accuracy of the model, polynomial and Bayesian optimization methods are used to calibrate the model. The density and scale control parameters are calibrated based on the boundary layer distortion prediction model. The disturbance parameters solved by the boundary layer linear stability theory are used as the motion parameters of the proxy model. The control parameters of the flow direction movement and dip distribution are calibrated based on the turbulent density fluctuation distribution law. Finally, a density proxy model reflecting the optical characteristics of high-speed flow field is obtained, and the fast simulation of aero-optical effects is completed through the calculation of light. The remainder of this paper is organized as follows: in Section 2, the design of the density proxy model in DP-AOQS is presented. The calibration method of the control parameters in the density proxy model is described in Section 3. In Section 4, the aero-optical effects in the visible light band are verified from the perspectives of density field distribution, OPD, and fuzzy star map. Finally, the conclusions are drawn in Section 5.

Design of Density Proxy Model
In this section, the density proxy model of the boundary layer flow field is established. This paper mainly studies the boundary layer flow field outside the optical window, and the thickness of the TBL under hypersonic conditions is about 5-20 mm. In addition, considering the size of the optical window, in order to fully simulate the transmission disturbance of the incident starlight by the TBL flow field outside the optical window fully, a cuboid space above the optical window was selected as the simulation domain with the size of  The density sphere vortex model was used to establish the turbulence density proxy model, in which the density sphere vortex model is used to simulate the large-scale The remainder of this paper is organized as follows: in Section 2, the design of the density proxy model in DP-AOQS is presented. The calibration method of the control parameters in the density proxy model is described in Section 3. In Section 4, the aerooptical effects in the visible light band are verified from the perspectives of density field distribution, OPD, and fuzzy star map. Finally, the conclusions are drawn in Section 5.

Design of Density Proxy Model
In this section, the density proxy model of the boundary layer flow field is established. This paper mainly studies the boundary layer flow field outside the optical window, and the thickness of the TBL under hypersonic conditions is about 5-20 mm. In addition, considering the size of the optical window, in order to fully simulate the transmission disturbance of the incident starlight by the TBL flow field outside the optical window fully, a cuboid space above the optical window was selected as the simulation domain with the size of L x × L y × L z = 200 mm × 30 mm × 100 mm, as shown in Figure 2. OXYZ is defined as the window coordinate system. The coordinate origin O is located at the centroid of the upper surface of the optical window. The remainder of this paper is organized as follows: in Section 2, the design of the density proxy model in DP-AOQS is presented. The calibration method of the control parameters in the density proxy model is described in Section 3. In Section 4, the aero-optical effects in the visible light band are verified from the perspectives of density field distribution, OPD, and fuzzy star map. Finally, the conclusions are drawn in Section 5.

Design of Density Proxy Model
In this section, the density proxy model of the boundary layer flow field is established. This paper mainly studies the boundary layer flow field outside the optical window, and the thickness of the TBL under hypersonic conditions is about 5-20 mm. In addition, considering the size of the optical window, in order to fully simulate the transmission disturbance of the incident starlight by the TBL flow field outside the optical window fully, a cuboid space above the optical window was selected as the simulation domain with the size of  The density sphere vortex model was used to establish the turbulence density proxy model, in which the density sphere vortex model is used to simulate the large-scale The density sphere vortex model was used to establish the turbulence density proxy model, in which the density sphere vortex model is used to simulate the large-scale structure in turbulence that is much larger than the wavelength of the incident light and has density spatial coherence.

Continuous Density Model of Ellipsoidal Vortex
The density ρ V (r,t) of the large-scale vortex structure is composed of multiple ellipsoidal vortex densities, which can be expressed as follows: where ρ i (r,t) is the density distribution of the i-th ellipsoidal vortex, which is different from the uniform ellipsoidal vortex in the literature [27]. In order to ensure the continuity of the density field, the internal density ρ i ( i r, i ψ, i h) of the ellipsoidal vortex model is designed as a layered distribution for time t: where i r, i ψ and i h are the cylindrical coordinates under the ellipsoidal vortex coordinate system o i x i,V y i,V z i,V , as shown in Figure 3. structure in turbulence that is much larger than the wavelength of the incident light and has density spatial coherence.

Continuous Density Model of Ellipsoidal Vortex
The density ( ) V t ρ r, of the large-scale vortex structure is composed of multiple ellipsoidal vortex densities, which can be expressed as follows: is the density distribution of the i -th ellipsoidal vortex, which is different from the uniform ellipsoidal vortex in the literature [27]. In order to ensure the continuity of the density field, the internal density ( , , ) of the ellipsoidal vortex model is designed as a layered distribution for time t : where , in ρ is the maximum spatial fluctuation in the density inside the ellipsoidal vortex, which can be expressed as follows: where ρ ∞ is the incoming flow density under the current conditions. ( ) 1 N − represents that the density change inside the ellipsoidal vortex may be higher or lower than the incoming flow density. Λ is the major axis length of the ellipsoidal vortex and k ρ is the control parameter of the ellipsoidal vortex density.

( )
, F M a y ∞ is the Mach number influence factor, which is given by: o i is the centroid of the i-th ellipsoidal vortex; axis x i,V points to the major axis of the ellipsoid and is located in plane OXY; axis o i y i,V points to the minor axis; and Λi h is the major axis length of the ellipsoidal surface at a distance i h from the plane o i x i,V y i,V . ρ in is the maximum spatial fluctuation in the density inside the ellipsoidal vortex, which can be expressed as follows: where ρ ∞ is the incoming flow density under the current conditions. (−1) N represents that the density change inside the ellipsoidal vortex may be higher or lower than the incoming flow density. Λ is the major axis length of the ellipsoidal vortex and k ρ is the control parameter of the ellipsoidal vortex density. F(Ma ∞ , y) is the Mach number influence factor, which is given by: where f (y/δ) = U y/δ /U ∞ is the average flow velocity distribution; U y/δ is the flow velocity at the height y/δ, which can be approximated by the logarithmic law c 1 ln(y/δ) + c 2 in a completely turbulent state; and U ∞ is the incoming flow velocity.

Location and Scale Model of Ellipsoidal Vortex
The position and scale distribution of the ellipsoidal vortex structure in the density proxy model are modeled by using the increasing function relationship. The position distribution of the ellipsoidal vortex is established as a Gaussian distribution model about the wall height: where y = y/δ and P y are the probability of arranging an ellipsoidal vortex at the wall height y. y 0 is the wall height with the maximum probability of placing the ellipsoidal vortex. Results from different studies have shown that structures at a height of approximately 0.8δ contribute most to optical distortions [14,28], so we took y 0 = 0.8δ. σ is the variance of Gaussian distribution, and we took it as a fixed value of 0.5δ and used the size of the gas ellipsoid to adjust the shape of the distorted wavefront. The normalization coefficient K satisfies: where k num is the total number of ellipsoidal vortices in the simulation domain, and the typical value is between hundreds and thousands [29], so k num = 5000 was taken in this paper. The size of the large-scale structure in the boundary layer increases at a certain height. The major axis length of the ellipsoidal vortex is described by a quadratic function as follows: Λ y = k len y − y min 2 + Λ min (7) where y min is the wall height of the smallest ellipsoidal vortex; namely, the height for the vertex of the quadratic function, taken as 0.05δ according to the literature [14]. Λ min is the minimum value of Λ y . The range of the minimum density correlation length in TBL is approximately 0.05δ∼0.1δ [14,28], so we took Λ min = 0.05δ. k len is taken as the control parameter of the gas-ellipsoidal scale. The density correlation length range of the turbulent boundary layer is about 0.15δ ≤ Λ(y) max ≤ 0.4δ [17,27]. The constraint range of k len can be obtained according to Equation (7).

Motion Model of Ellipsoidal Vortex
There are disturbance waves with different wave numbers and phase velocities in the supersonic boundary layer [30]. The motion of the ellipsoidal vortex model is modeled as disturbance waves, in which the centroid (x i , y i , z i ) of the ellipsoidal vortex moves in the form of the superposition of the flow direction traveling waves: whereû k is the traveling wave disturbance characteristic function; α k is the k-th order disturbance flow beam; β k is the k-th order disturbance spread beam; ω k is the k-th order disturbance frequency; ϕ k is the k-th order mode phase; and j is an imaginary unit. For the time mode, α k is a real number and ω k is a complex number; that is, ω k = ω k,r + ω k,j j. c.c is a complex conjugate.
In the hypersonic boundary layer, the first mode has an important influence on the transition, and the second mode disturbance becomes the most unstable disturbance wave after the incoming Mach number exceeds four. Therefore, the ellipsoidal vortex center motion is set as the superposition of the first mode and the second mode; that is, k = 2. The first mode flow direction beam α 1 is 0.10892184, the spread direction beam β 1 is ±0.313, the second mode flow direction beam α 2 is 1.44368594, and β 2 is ±1.5. According to the linear stability theory, the disturbance characteristic functionsû k and ω k in the control parameters of the ellipsoidal vortex motion can be obtained by solving the Orr-Sommerfeld small-disturbance equation. Considering the directional dependence of the large-scale turbulent structure, the expected inclination of the ellipsoidal vortex was set to 24.17 • , and the oblateness was set to 0.447.

Calibration of Control Parameters in Density Proxy Model
In this section, the control parameters of the density proxy model are calibrated. k ρ , k len was used to control the position distribution of the ellipsoidal vortex and the spatial fluctuation of the density in the density proxy model, which can be calibrated by the methods of polynomial and Bayesian optimization.

Density and Scale Control Parameter Constraints
The density control parameter k ρ affects the density spatial fluctuation of the ellipsoidal vortex structure, and the scale control parameter k len affects the size distribution of the ellipsoidal vortex, both of which are closely related to the incoming flow state. The density correlation length range of the turbulent boundary layer is about 0.15δ∼0.4δ [17,27]. The constraint range of k len can be obtained according to Equation (7). To determine the constraint range of k ρ , the Strehl Ratio (SR) of the density proxy model was calculated in a steady state for different control parameters under the condition of a flight altitude of 10 km and Mach number of 5 Ma, as shown in Figure 4.   sulting in the excessive spatial pulsation of the turbulence proxy model. Therefore, the upper limit b k ρ and lower limit a k ρ are used for the density control parameter k ρ , meeting the following conditions:  It can be seen from Figure 4 that SR s gradually decreases from 1 to 0 with the increase in the control parameters k ρ and k len . When k ρ ≤ 0.001, the change in SR s with the increase in k len can be ignored, and the turbulence density proxy model cannot reflect the spatial fluctuation of turbulence. When k ρ ≥ 100, SR s decreases rapidly with the increase in k len and the decreasing speed continues to increase with the increase in k ρ , resulting in the excessive spatial pulsation of the turbulence proxy model. Therefore, the upper limit k b ρ and lower limit k a ρ are used for the density control parameter k ρ , meeting the following conditions: where SR a and SR b are the minimum and maximum values of SR under the working conditions to be calibrated, respectively.

Calibration of Density and Scale Control Parameters
Under the given inflow conditions, the mean square value of OPD of the TBL can be described by the boundary layer distortion prediction model [25]: where B is a constant, taken as 0.2, F(Ma ∞ ) is an empirical function of the incoming Mach number, and C f is the wall friction coefficient of the TBL of the flat plate, which is related to the Reynolds number Re x . The density fluctuation and scale structure in turbulence are affected by the Mach number and Reynolds number of the incoming flow. The density control parameters affect the density fluctuation of the ellipsoidal vortex structure, and the scale control parameters affect the size distribution of the ellipsoidal vortex, both of which determine the optical characteristics of the density proxy model. Therefore, the control parameters k ρ and k len can be calibrated as a function of Ma ∞ , Re x , and flight altitude H, which can be expressed as follows: The working condition is expressed as y = Ma ∞ Re x H T .g len (Ma ∞ , Re x , H) and g ρ (Ma ∞ , Re x , H) are expressed as g(y). With the boundary layer distortion prediction model as the constraint condition, the optimal function mapping g * is obtained; that is, a variational problem is solved: where OPD p rms (y) is the root mean square of the OPD of the boundary layer distortion prediction model under the given working condition y, and OPD s rms (g(y)) is the root mean square of the OPD of the density proxy model. Because the gradient information of the control parameter output by the density proxy model cannot be obtained, the traditional variational method cannot be used to solve the variational problem in Equation (12). We transformed it into a parameter optimization problem, which is solved by a gradient-free parameter optimization algorithm.
SR in a density proxy model and distortion prediction model can be described as the form of multiplication or the addition of two univariate continuous functions with variable separation. Considering the existence of the multiplication form of the Mach number and Reynolds number in the distortion prediction model, we conducted research in the form of multiplication: where x 1 and x 2 in the density proxy model are k len and k ρ , respectively, and h 1 (x 1 ) and h 2 (x 2 ) are the monotonic increasing functions. x 1 and x 2 in the distortion prediction model are Ma ∞ and Re x , respectively, h 2 (x 2 ) is the monotonic decreasing function, and the monotonicity of h 1 (x 1 ) is uncertain. Therefore, the control parameters k len and k ρ can be calibrated as the multiplication of two continuous functions with one variable separated: The univariate continuous functions g len,i and g ρ,i can be fitted by the power function polynomial. The polynomial fitting function is set as g(x; a, N), where N is the highest order term of the power function and a = a 0 a 1 · · · a N T is the coefficient vector of the power function.
g(x; π, N) = ∑ N k=0 a k x k The coefficient a 1 , a 2 , a 3 , a 4 of the power function polynomial and the highest degree term N 1 , N 2 , N 3 , N 4 of the power function polynomial are parameters to be calibrated. For the convenience of writing, A and N are represented as A = a 1 a 2 a 3 a 4 and T , respectively. Then, Equation (11) can be written as: The surface shape of SR p distribution with Ma ∞ and Re x in the distortion prediction model changes obviously with the height. When fitting the functional relationship between k ρ , k len and {Ma ∞ , Re x } by polynomial, it is necessary to calibrate the power function coefficient matrix A and the highest-order term N by height H. In order to determine the range of coefficient A to be optimized, the input variable Ma ∞ , Re x , and the output variable k ρ , k len are normalized, and the value range is limited to the value interval [0, 1]. Ma ∞ and log(Re x ) are calibrated according to the range of working conditions. Considering the scope of application of the distortion prediction model, the value range of Ma ∞ is limited to [1.5, 5]; Re x is limited to 10 6 , 10 8 ; k len is limited to k a len /γ R , k b len /γ R ; and k ρ is limited to k a ρ /γ R , k b ρ /γ R . γ R is a trust coefficient which is used to change the range of the control parameters to be calibrated. Additionally, γ R is considered as a hyper-parameter, and the value range is set to [0. 5,2]. When γ R < 1, the calibration range of the control parameters is reduced, and when γ R > 1, the calibration range of the control parameters becomes larger.
After the normalization of the input variables Ma ∞ , log(Re x ) and output variables k len , k ρ , the value of the higher-order term a ij of the power function polynomial coefficient is naturally limited to [−1, 1], and the range of the zero-order power function coefficient A is [0, 1]. Thus, coefficient matrix A satisfies the following equation: The minimum value of the highest power term of the constrained power function is 3 and the maximum value is 20. The parameters A and N are calibrated by iterative optimization. For the fixed hyper-parameters N and γ R at the given flight altitude H 0 , the parameter A is optimized and solved. With the boundary layer distortion prediction model as the constraint, the optimal solution A * of the parameter A to be calibrated is an optimization problem: where OPD s rms (Ma ∞ , Re x ; A) is the root mean square of the simulated OPD of the density proxy model at flight altitude H 0 . For a given flight altitude H 0 , the parameter N and trust coefficient γ R can be regarded as the hyper-parameters for optimization, and the loss function is defined as: where n iter is the total number of iterations in the optimization, a P log(n iter ) is the penalty function of the total number of iterations, and a P is the penalty coefficient, taken as 0.001. When the number of iterations is too large, it means that the iterative optimization effect of the polynomial parameters is poor, and the polynomial form is not suitable as a fitting function. For the above hyper-parameter optimization problem, the Bayesian optimization algorithm is used to solve it: (20) where x is the parameter to be optimized, and its value range is X. f (x) is the objective function, and its gradient information is unknown. The basic idea of Bayesian optimization is to estimate the posterior distribution of the objective function based on the Bayesian principle. It is usually necessary to make assumptions about the model of the posterior distribution, estimate the posterior distribution on the basis of the observed data, and select new sampling points by maximizing the effect function to achieve the posterior distribution approaching the true distribution [31].
For H 0 = 15 km, the solution results of the hyper-parameters based on Bayesian optimization are shown in Figure 5, where (a) and (b) are the iterative scatter distribution of the hyper-parameter N. The red point is the optimal value, and the blue point is the initial iteration point. The size of the point increases with the number of iterations. Figure 5c shows the iteration trajectory of the trust coefficient γ R , which reaches the optimal value after 30 iterations. The change curve of the loss function Loss is shown in Figure 5d. The prismatic sign is the loss function observed at the sampling point in the current iteration. The circle sign is the minimum value of the loss function in the current iteration. When the number of iterations is 12, the minimum value of the loss function is close to convergence. When the minimum value of the loss function reaches the tolerance range in 30 iterations, the Bayesian optimization solution is completed.  SR s of the turbulence density proxy model and SR p of the distortion prediction model after calibration are shown in Figure 6a,b. After polynomial calibration, the fitting error |SR s − SR p | of SR is within 2 × 10 −3 .

Simulation and Analysis
In this section, with a star sensor as the object, the simulation performance of DP-AOQS is verified from the characteristics of the density field distribution, OPD, and fuzzy star map of the density proxy model.

Verification of Optical Characteristics of Density Proxy Model
For the simulation conditions of . With the passage of time, the movement of the vortex structure in the direction of flow made the distribution of OPD fluctuate in both the time and space dimensions, which is consistent with the pulsation characteristics of the turbulent boundary layer [32]. It only took 200 s for DP-AOQS to obtain a set of results (simulation environment: Matlab R2021b, Intel i9-9900X CPU, 128 G RAM). A high-precision flow field structure calculated by CFD required 10-20 h (simulation environment: a high-performance server with AMD EPYC ROME 7H12, 128-cores CPU, 256 G RAM, the number of grids was more than 6 million, and the generation of large eddy simulation (LES) was more than 20000).

Simulation and Analysis
In this section, with a star sensor as the object, the simulation performance of DP-AOQS is verified from the characteristics of the density field distribution, OPD, and fuzzy star map of the density proxy model.

Verification of Optical Characteristics of Density Proxy Model
For the simulation conditions of H = 15 km, Ma ∞ = 3, Re x = 10 6 , the density fluctuation, and OPD of the density proxy model in the XY plane after being calibrated in Section 3 are shown in Figure 7. It can be seen from Figure 7 that the OPD fluctuated greatly with the flow direction distance because there were ellipsoidal vortices with uneven density fluctuations in the flow direction. For example, at t = 5 ms, there were wave peaks in the OPD near x = −0.08 m, and there were many ellipsoidal vortices (near y = 0.02 m − 0.03 m) with positive density fluctuations near x = −0.08 m. With the passage of time, the movement of the vortex structure in the direction of flow made the distribution of OPD fluctuate in both the time and space dimensions, which is consistent with the pulsation characteristics of the turbulent boundary layer [32]. It only took 200 s for DP-AOQS to obtain a set of results (simulation environment: Matlab R2021b, Intel i9-9900X CPU, 128 G RAM). A high-precision flow field structure calculated by CFD required 10-20 h (simulation environment: a high-performance server with AMD EPYC ROME 7H12, 128-cores CPU, 256 G RAM, the number of grids was more than 6 million, and the generation of large eddy simulation (LES) was more than 20,000).
In order to compare with the OPD calculated by the CFD, the OPD and projection surface of the density proxy model were obtained under the condition H = 20 km, Ma ∞ = 3.8, as shown in Figure 8a. The OPD of the flow field near the concave window under the same simulation conditions by the CFD is shown in Figure 8b [19]. It can be seen that there were also many peaks and troughs in the optical path difference surface. Figure 8 demonstrates that DP-AOQS can reflect the variation trend of the optical properties of the actual turbulent boundary layer flow field.
The flow direction OPD caused by the density fluctuation of the turbulent structure has certain amplitude frequency characteristics. Through experiments and similarity analysis, the amplitude spectrum model of the flow direction OPD was determined as follows [25]: where St δx is the normalized frequency of flow direction, satisfying St δx = δ f s /U ∞ . f s is the sampling frequency, ρ SL is the reference density of the sea level, ρ ∞ and U ∞ are the incoming flow density and velocity, respectively, andθ peak is the peak value of the amplitude. To eliminate the difference caused by working conditions, we normalized Equation (21) to obtain the flow direction normalized spectrum OPD x (St δx ), as shown in the following equation: In order to compare with the OPD calculated by the CFD, the OPD and projection surface of the density proxy model were obtained under the condition 20km, 3.8 H M a ∞ = = , as shown in Figure 8a. The OPD of the flow field near the concave window under the same simulation conditions by the CFD is shown in Figure 8b [19]. It can be seen that there were also many peaks and troughs in the optical path difference surface. Figure 8 demonstrates that DP-AOQS can reflect the variation trend of the optical The flow direction OPD caused by the density fluctuation of the turbulent structure has certain amplitude frequency characteristics. Through experiments and similarity analysis, the amplitude spectrum model of the flow direction OPD was determined as follows [25]: The amplitude spectrum model, density proxy model, and CFD flow field in the literature [19] were used to calculate the OPD spectrum distribution, as shown in Figure 9. It can be seen from Figure 9 that the density proxy model in this paper had a downward trend as the flow direction frequency increased, which is consistent with the trend of the amplitude spectrum model and the OPD spectrum obtained from the CFD calculation. However, the amplitude of the OPD spectrum in this paper decreased rapidly in the interval 10 10 x St δ = − , which is caused by spatial discrete sampling, and amplitude mutation will occur at uncovered frequency points. The amplitude spectrum model, density proxy model, and CFD flow field in the literature [19] were used to calculate the OPD spectrum distribution, as shown in Figure 9. It can be seen from Figure 9 that the density proxy model in this paper had a downward trend as the flow direction frequency increased, which is consistent with the trend of the amplitude spectrum model and the OPD spectrum obtained from the CFD calculation. However, the amplitude of the OPD spectrum in this paper decreased rapidly in the interval St δx = 10 −1 − 10 0 because the flow direction motion model only used the first mode and the second mode for superposition, while there are many higher-order mode disturbance waves in the boundary layer. It is normal for the OPD spectrum to show buffeting in the interval St δx = 10 0 − 10 1 , which is caused by spatial discrete sampling, and amplitude mutation will occur at uncovered frequency points.

Monte Carlo Simulation under Different Working Conditions
The OPD of the density proxy model at 10 Figure 10, where there were 100 sampling points in total.

Monte Carlo Simulation under Different Working Conditions
The OPD of the density proxy model at 10 km and 15 km was verified by the Monte Carlo method, where the condition was Ma ∞ ∈ [2,5], Re x ∈ 1 × 10 6 , 1 ×  From the relative error distribution of OPD in Figure 10, the fitting error of the density proxy model was within 10% at 15 km, while the flight altitude was 10 km, and the fitting error was within 12%, satisfying the accuracy requirements of the simulation. When the altitude decreases, the fitting error increases. This is because when the height decreases, the surface of SR will bend more violently with {Ma ∞ , Re x }, and the error of the control parameters at the uncalibrated height increases. In general, the calibration of the control parameters of the density proxy model at a given altitude within 1 km could satisfy the accuracy requirements of OPD, and the calibration accuracy was higher than the steadystate simulation method proposed in the literature [27].

Simulation of Distorted Star Maps
After the starlight passes through the high-speed flow field, the star map on the image plane of the star sensor will degenerate. Because the distance between the star and the aircraft observed by the star sensor is far greater than the distance between the Earth and the sun, the observed star is regarded as a point light source at infinity. The star point signal on the CCD image plane is concentrated in a small circular area, where the energy distribution can be assumed to be the Gaussian distribution. The light field distribution on the image plane is the gray distribution of the diffraction image of a single star, and the display gray level of the corresponding pixel is obtained after gray level aggregation.
Under the condition of Ma ∞ = 3, Re x = 1 × 10 6 , the simulation results of the two distorted stars at 15 km and 5 km were as shown in Figure 11. In the simulation, the focal length of the star sensor was f = 61.88 mm, the pupil size was 50 mm, the pixel size was 15µm, the field angle was FOV = 20 • , and the number of pixels on the image plane was 1456 × 1456. the display gray level of the corresponding pixel is obtained after gray level aggregation.
Under the condition of 6 3, 1 10 x Ma Re ∞ = = × , the simulation results of the two distorted stars at 15 km and 5 km were as shown in Figure 11. In the simulation, the focal length of the star sensor was 61.88mm f = , the pupil size was 50 mm, the pixel size was 15μm, the field angle was 20 FOV = ο , and the number of pixels on the image plane was  The gray value at the center of the distorted star map is listed in Table 1. When the flight altitude was 15 km, the SR p of the distortion prediction model was 0.8613, and the SR s of the turbulence density proxy model was 0.8650. At this point, compared with the standard star map, the gray distribution of star A and star B in the distorted star map was still the center point, and the energy was mainly concentrated in the range of one pixel in each direction from the center point. The gray value of the center point in the distorted star map decreased by 5-9 gray values. When the flight altitude was 5 km, the SR p of the distortion prediction model was 0.0621, and the SR s of the turbulence density proxy model was 0.0645. Compared with the standard star map, the brightest point is still the center point of stars A and B in the distorted star map, but the energy distribution has spread in all directions, spreading 1-2 pixels. The gray value of the center point in the distorted star maps is greatly reduced, and the gray value of star A is the most obvious, falling 79 gray values. The distorted star maps obtained by DP-AOQS can be used for a reference in the design stage of celestial navigation in hypersonic vehicles. In addition, in order to highlight the advantages of this study, the distorted star map results of this paper are compared with the previous research methods. The star map obtained from the previous research [27] is shown in Figure 12. pixels. The gray value of the center point in the distorted star maps is greatly reduced, and the gray value of star A is the most obvious, falling 79 gray values. The distorted star maps obtained by DP-AOQS can be used for a reference in the design stage of celestial navigation in hypersonic vehicles. In addition, in order to highlight the advantages of this study, the distorted star map results of this paper are compared with the previous research methods. The star map obtained from the previous research [27] is shown in Figure  12. Figure 12. The simulation results of distorted star maps with the previous study [27].
By comparing the star map results in Figure 12 with those in Figure 11c,d, it can be found that the star map results are in line with objective laws, and the diffusion effect of the star map is more obvious because of the density change in the interior of the ellipsoid vortex and the motion model of the ellipsoid vortex considered in this paper.

Conclusions
In this paper, a quick simulation method of the aero-optical effects based on the density proxy model (DP-AOQS) is proposed. The density proxy model, considering the optical characteristics of the turbulent steady state and motion state, is designed to replace the density field in CFD numerical simulation. Based on the aero-optical laws obtained Figure 12. The simulation results of distorted star maps with the previous study [27].
By comparing the star map results in Figure 12 with those in Figure 11c,d, it can be found that the star map results are in line with objective laws, and the diffusion effect of the star map is more obvious because of the density change in the interior of the ellipsoid vortex and the motion model of the ellipsoid vortex considered in this paper.

Conclusions
In this paper, a quick simulation method of the aero-optical effects based on the density proxy model (DP-AOQS) is proposed. The density proxy model, considering the optical characteristics of the turbulent steady state and motion state, is designed to replace the density field in CFD numerical simulation. Based on the aero-optical laws obtained from the experiments and simulation analysis, the density proxy model is calibrated by combining polynomial fitting and Bayesian optimization so that the proxy model can reflect the optical characteristics in the actual high-speed flow field. The performance of the aero-optical effects in the visible light band is verified from the perspectives of the density field distribution, OPD, and fuzzy star map. Compared with the traditional CFD simulation method, DP-AOQS can quickly provide the distortion results of aero-optical effects under different flight conditions. It provides a powerful tool for aero-optical effects analysis in the design stage of celestial navigation in hypersonic vehicles.   γ R Trust coefficient which is used to change the range of control parameters.

Loss
Loss function. n iter Total number of iterations in the optimization. a P log(n iter ) Penalty function of the total number of iterations. a P Penalty coefficient.